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Abstract. This work evaluates the interaction of a simulated spectrum of 
convectively generated gravity waves with realistic middle atmosphere mean winds. 
The wave spectrum is derived from the nonlinear convection model described by 
Alexander et al [1995] that simulated a two-dimensional midlatitude squall line. 
This spectrum becomes input to a linear ray tracing model for evaluation of wave 
propagation as a function of height through climatological background wind and 
buoyancy frequency profiles. The energy defined by the spectrum as a function 
of wavenumber and frequency is distributed spatially and temporally into wave 
packets for the purpose of estimating wave amplitudes at the lower boundary of the 
ray tracing model. A wavelet analysis provides an estimate of these wave packet 
widths in space and time. Without this redistribution of energies into wave packets 
the Fourier analysis alone inaccurately assumes the energy is evenly distributed 
throughout the storm model domain. The growth with height of wave amplitudes 
is derived from wave action flux conservation coupled to a convective instability 
saturation condition. Mean flow accelerations and wave energy dissipation profiles 
are derived from this analysis and compared to parameterized estimates of gravity 
wave forcing, providing a measure of the importance of the storm source to global 
gravity wave forcing. The results suggest that a single large convective storm 
system like the simulated squall line could provide a significant fraction of the 
zonal mean gravity wave forcing at some levels, particularly in the mesosphere. 
The vertical distributions of mean flow acceleration and energy dissipation do not 
much resemble the parameterized profiles in form because of the peculiarities of 
the spectral properties of the waves from the storm source. The ray tracing model 
developed herein provides a tool for examining the role of convectively generated 
waves in middle atmosphere physics. 


1. Introduction 

Gravity waves transport energy and momentum from 
the troposphere to the middle atmosphere where, it is 
widely recognized, they can have a profound effect on 
the general circulation patterns, temperature structure, 
and spatial distributions of mixing ratios of the atmo- 
spheric gases. The importance of wave drag and diffu- 
sion in the middle atmosphere was clearly demonstrated 
in zonal mean model studies in the 1980s [e.g. Holton, 
1982, 1983; Dunkerton, 1982; Garcia and Solomon , 
1985] that utilized the gravity wave parameterization 
developed by Lindzen [1981]. Lindzen’s parameteriza- 
tion required assumptions about the phase speeds and 
source distributions of gravity waves which are to date 
still not well characterized. Holton’s [1982] work, as well 
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as that of Matsuno [1982], established the importance of 
high phase speed waves to explain the observed mean 
zonal wind structure and thereby stressed the impor- 
tance of wave sources other than flow over topography. 

These wave-driven processes are also important in 
three-dimensional global circulation models where pa- 
rameterization of gravity wave effects is complicated by 
the models’ sensitivity to the geographical and tempo- 
ral distributions of wave sources. Planetary scale waves 
can be resolved explicitly in these models, however the 
effects of smaller-scale waves (of the order of 100 km 
and less) will likely be treated only via parameterization 
for some time to come. Orographically excited waves 
have been successfully parameterized in such models 
and have been shown to affect circulation in the tropo- 
sphere as well as the middle atmosphere [Palmer et a/., 
1986; McFarlane , 1987; Bacmeister , 1995]. Other wave 
sources have been more difficult to characterize. Specif- 
ically, waves excited by convective activity are likely 
very important in the tropics and southern hemisphere 
where there are few orographic wave sources, and con- 
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vection may be a source of the high phase speed waves 
known to be important in the mesosphere. 

Recent modeling efforts reported by Fovell et al. 
[1992], Holton and Durran [1993], and Alexander et al. 
[1995] have described properties of vertically propagat- 
ing waves generated by deep convection in the form of 
a two-dimensional squall line. In the work of Alexander 
et al [1995; hereafter referred to as AHD] the spec- 
tral properties of the waves in the stratosphere in their 
simulation were characterized and compared to observa- 
tions. They noted a strong response at high frequencies 
(corresponding to periods of 8 min to 1 hour) and at 
long vertical wavelengths (6-10 km), which show sim- 
ilarity to observations of stratospheric motions above 
convective sources [Larsen et a/., 1982; Sato , 1993], and 
which may be characteristic of waves associated with 
deep convection. Additional modeling studies and ob- 
servations will eventually clarify the generality of these 
spectral characteristics. 

In this work, the interaction of this spectrum of con- 
vectively generated waves with realistic middle atmo- 
sphere winds is examined via a linear ray tracing tech- 
nique which uses conservation of wave action flux to 
predict wave amplitudes as a function of height. Wave 
interactions with the mean flow are included as satu- 
ration effects when amplitudes exceed convective insta- 
bility limits and via the filtering effects of critical level 
absorption and wave reflection. This method is chosen 
as a means of estimating the importance of the contri- 
bution of waves from such a convective source to the 
estimated global gravity wave forcing. It is unrealis- 
tic to extend the domain of the nonlinear simulation of 
AHD to the altitudes and horizontal distances required 
to make this estimate. Durran [1995] also highlights 
the difficulty in using traditional diagnostic methods for 
evaluating gravity wave dissipation effects in a domain 
of limited size. The approach here is to take the two- 
dimensional power spectrum derived from the nonlin- 
ear model (AHD) as input to a linear wave propagation 
calculation. The power in this spectrum is distributed 
into packets of finite width in horizontal distance, z, 
and time, t, and the packet widths are estimated from a 
wavelet analysis of the nonlinear model results. The lin- 
ear and nonlinear results cure tested for consistency be- 
low 32-km altitude where the models overlap. The ray 
tracing model is in many ways simpler than the three- 
dimensional models described by Eckermann [1992] and 
Marks and Eckermann [1995] that were designed to 
study global propagation characteristics over the full 
range of possible gravity wave frequencies. For the spec- 
trum of convectively generated waves considered here, a 
number of simplifying assumptions are possible for the 
purpose of estimating the mean flow forcing. 

The results of the linear wave propagation include 
estimates of the mean flow acceleration and energy dis- 
sipation rates associated with the input wave spectrum 
and specified climatological mean wind profiles. Com- 
parison to the spectral gravity wave parameterization 
of Fritts and Lu [1993] for the same mean wind profiles 
provides an estimate of the importance of the single 


storm source to global gravity wave forcing and also 
gives insight into how the peculiar spectral character- 
istics associated with the convectively generated waves 
affect the profile of wave drag and dissipation. The 
method developed here can be used to compare future 
simulations to the midlatitude case of AHD and also 
provides an avenue for testing simpler parameteriza- 
tions of convectively generated waves against this more 
complete spectral description. 

The following section briefly reviews the convection 
simulation and determination of the two-dimensional 
power spectrum described by AHD. The method of con- 
verting power spectral density to wave amplitude is also 
derived, including a wavelet analysis of the nonlinear 
model results to estimate wave energy packet dimen- 
sions in space and time. This then establishes lower 
boundary conditions for the linear ray tracing analysis 
described in section 3, which includes determinations 
of wave amplitudes as a function of height as well as 
effects on the mean state. Expressions for mean flow 
acceleration and rate of wave energy dissipation are de- 
rived. Section 4 is devoted to checking the linear prop- 
agation model against the full nonlinear model results 
between 13 and 32 km where the two models overlap, 
testing many of the simplifying assumptions in the lin- 
ear model. In section 5 the wave spectrum interaction 
with realistic mean wind and buoyancy frequency pro- 
files is examined and compared to the Fritts and Lu 
[1993] parameterized results for the same background 
state. A concluding summary follows in section 6. 

2. Determination of Gravity Wave 
Amplitudes and Propagation 
Characteristics 

2.1. The Convection Simulation 

The convectively generated gravity wave spectrum in 
this analysis is derived from the results of a numerical 
midlatitude squall line simulation previously described 
by AHD [Alexander et al ., 1995], Some of the features 
of the model will be briefly described here. For more 
details the reader is also referred to previous applica- 
tions of the squall line model described by Fovell et al. 
[1992] and Holton and Durran [1993]. 

The convection simulation from which the gravity 
wave spectrum is derived is a two-dimensional, nonlin- 
ear, compressible, nonhydrostatic squall line model. It 
resolves a deep stratosphere layer, from the tropopause 
at ~12 km up to 32 km altitude. The full model domain 
is 840 km in the horizontal and 32 km in the vertical and 
includes wave permeable boundaries at the sides and 
top. The model reference frame translates eastward at 
16ms" 1 to track the motion of the storm, keeping it in 
the center of the domain. Winds in the stratosphere in 
this model are constant in height at 16 m s -1 , the same 
as the reference frame translation speed, and so the 
stratospheric waves produced by the storm are viewed 
in their “intrinsic” frame of reference (i.e., the observed 
wave frequencies are the intrinsic frequencies). A rich 
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Figure 1. Contours of vertical velocity in the stratospheric portion of the nonlinear convection 
model. Contours are plotted at 0.5 m s' 1 intervals. Dotted contours represent negative values, 
(a) w(x.z) at t = 5 hours. Horizontal lines represent surfaces of constant potential temperature 
at 25 K intervals. The storm center, which is the main source region for these waves, lies in the 
troposphere below the figure and is located at x ~ 420 km. (b) w(x, t) at z = 30 km. The slopes 
of surfaces of constant phase indicate the intrinsic phase speeds of the wave motions. 


spectrum of gravity waves is generated in the strato- 
sphere, as can be seen in Figure 1. Figure la shows a 
single time frame of vertical velocity contours and po- 
tential temperature surfaces above 13 km. In Figure 
lb, contours of w(x,t) at z = 30 km are shown. The 
slope of the phase surfaces in Figure lb is an indication 
of the separation east and west of storm center (at x ~ 
420 km) of eastward and westward propagating waves, 
respectively, and points to the storm center region as 
the location of the primary source of the wave energy. 
A strong preference for forcing of westward propagat- 
ing waves is observed in this and other similar simu- 
lations due to a westward tilt with height of the main 
tropospheric updraft and the westward propagation of 
convection cells in the troposphere. These features have 
been described in earlier work with this model [Fovell et 
al, 1992; AHD], and have been observed in squall lines 
in nature [e.g. Home, 1993]. A spectral analysis of the 
stratospheric waves was described by AHD and revealed 
some distinctive spectral signatures which may be char- 
acteristic of waves generated by deep convection. AHD 
describe wave forcing mechanisms and their signatures 
in the spectral response and also review similarities to 
observations of wave motions above convective sources. 


2.2. Fourier Spectral Analysis 

The two-dimensional power spectrum in horizontal 
wavenumber and intrinsic frequency P(k 9 u>) computed 
from the simulation results of AHD will be used as lower 
boundary input to the linear propagation analysis to 
follow. The spectrum shown in Figure 2 separates 
power into eastward and westward phase with positive 
and negative frequencies, respectively. The power in 
this spectrum defines wave energies at an altitude of 
13 km. The superimposed white contour surrounds re- 
gions of the spectrum with power greater than or equal 
to 5 x 10 4 (m/s) 2 (cycle/m)" 1 (cycle/s)"” 1 . Regions out- 
side this contour with power less than this threshold will 
be excluded from the following analysis. The contour 
includes over 93% of the total energy in the spectrum, 
and spectral points outside this region are very likely 
heavily contaminated by spectral bias from regions of 
high power. Each pixel inside the white contour will 
be treated as a monchomatic wave packet. The wave 
amplitude can be derived from (1) the power spectral 
density at that point, together with (2) the definition 
of the wave packet width in time and space, and (3) an 
estimate of power aliasing to other regions of the spec- 
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where W ku is the discrete Fourier transform of w(x, t ) 
at wavenumber A: and frequency w. The power at (jfe,u>) 
describes the mean square amplitude of the wave with 
those characteristics averaged over the ( x,t ) domain, 
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Figure 2. Two-dimensional power spectrum of the ver- 
tical velocity field as a function of horizontal wavenum- 
ber and intrinsic frequency. The spectrum is derived 
from the two-dimensional Fourier transform of the field 
w{x y t) like that shown in Figure lb. Spectra at each 
altitude in the model have been averaged by using an 
exp [-(z - 13 km)/ff] weighting factor such that the 
spectral power represents amplitudes at 13 km near the 
tropopause. 

trum. Parameter (2) will be estimated from a wavelet 
analysis of the vertical velocity field described in the 
next section. 

The power spectrum in Figure 2 resolves wavenum- 
bers and frequencies, 
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with As = 1.5 km, N x = 256, At = 2 min and N t = 
128. So pixel dimensions are Ak = 2.6 x 10 ” 6 cycle 
m 1 and A& ~ 6.5 x 10” 5 cycle s” 1 . Wave energy with 
wavenumbers and frequencies smaller than the mini- 
mum nonzero values ( k < Ak and u> < Aw) is also 
omitted from the linear propagation analysis but repre- 
sents less than 0.2% of the total energy in the spectrum. 

The power spectral density in Figure 2 can be re- 
lated to the vertical velocity amplitude (in m s” 1 ) if it 
is assumed that a given pixel represents a wave mode 
with the frequency and wavenumber associated with 
that pixel. The two-dimensional power spectral density 
Pkw is given by 


This amplitude would be exactly the wave amplitude 
if the assumptions (inherent in Fourier analysis) of sta- 
tionarity throughout the (®,£) domain and periodicity 
on the x and t intervals were both satisfied. On the con- 
trary, both of these assumptions are violated, as can be 
seen in Figure 1. There is significant spatial and tempo- 
ral variation in the spectral properties. The waves can 
be thought of as concentrated in wave packets in both 
space and time; i.e., the areal extent of a given mode 
over the (®,£) domain is limited in both dimensions. 
These effects make the amplitude predicted from (2) 
smaller than the true amplitude at any point in (x,t). 
The power Pk w in equation (2) must be adjusted by two 
multiplicative factors: 

1. The first factor is equal to the areal extent of 
the whole domain divided by the energy-weighted areal 
extent of the wave packet. This factor arises because the 
spectrum gives a measure of the mean square amplitude 
in the domain, not the amplitude in the wave packet. 

2. The wave packet envelope can be thought of as 
a taper function which creates bias in the spectrum, 
reducing the value of the power at the peaks in the 
spectrum and spreading that power over a broad range 
of wavelengths. For a rectangular wave packet the bias 
would be described by the Fejer kernel [Percival and 
Walden } 1993, Section 6.4]. The power at the central 
peak of the kernel is proportional to the length of the 
window, so this second factor is also proportional to 
the ratio of the domain area to the envelope area, just 
as in factor 1. The proportionality constant will vary 
between 1 and 2 depending on the shape of the wave 
packet envelope. A value of 1 corresponding to a rect- 
angular shape is assumed for simplicity. 

Thus if the wave packets were rectangular in shape, 
and covered an area n x n t within the entire domain (area 
NxN t ), then the wave amplitude within that packet be- 
comes 


Aku, = \lPkti> f 2 A&Aw 

L \ n x n t ) 


1/2 


( 3 ) 


with one factor of ( N x N t /n x n t ) from each of factors 
1 and 2. This equation will be used to estimate am- 
plitudes at the lower boundary for the purpose of pre- 
dicting breaking levels for each mode in the ray tracing 
analysis. Parseval’s theorem, however, demands that 
the total power in the spectrum equal the power in the 
original signal integrated over the (x, t) domain, 
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Figure 3. An example of the wavelet basis functions 
centered in the analysis domain (x, t). Peaks and valleys 
are plotted with solid and dotted contours, respectively. 
This is the 96-km wavelength and 32-min period case. 
The panels along the top and side show cross sections 
in x and t through the center. 

— Ax At 2h*.*)i s - ( 4 ) 

k,u *>< 

Therefore the wave packet amplitude adjustments in 
(3) must be removed prior to evaluating any net ef- 
fects of the waves on the mean flow to conserve en- 
ergy. Wave packet dimensions will be estimated via the 
wavelet analysis described below. 

2.3. Wavelet Analysis 

Wavelet analysis is a spectral analysis technique that 
retains information about the spatial/temporal location 
of variations in the spectral power. Without this infor- 
mation the amplitudes derived solely from the Fourier 


analysis cannot be used to predict realistic breaking lev- 
els. For an energy-conserving orthogonal wavelet set, 
the spatial/temporal information is gained at the ex- 
pense of spectral resolution: Given a 128-point time 
array, the number of resolved frequencies in the wavelet 
spectrum is only log2(128) = 7, compared to the 64 
nonzero frequencies obtained from a Fourier analysis. 
Therefore wavelet analysis cannot replace the Fourier 
analysis for our purposes but can provide a rough esti- 
mate of the spatial/temporal extent of a given wavenum- 
ber/frequency signal and thus an estimate of an effective 
wave packet width. 

The 512x128 array of w(x,t) with 1.5-km and 2- 
min resolution, as shown in Figure lb, is examined 
via wavelet analysis. A two-dimensional, orthogonal 
wavelet transform is applied and an energy spectrum 
computed as the square of the resulting wavelet coef- 
ficients. Figure 3 shows an example of the wavelet 
basis functions employed here. The wavelet transform 
utilizes the set of Daubechies wavelet filters with 20 
coefficients summarized in the Numerical Recipes sub- 
routines “pwtset” and “pwt” [Press et a/., 1993, sec. 
13.10]. The set of basis functions employed in the anal- 
ysis consists of translations and dilations/contractions 
of the function shown in Figure 3. The wavelengths 
and periods resolved in the analysis are summarized in 
Table 1, as well as the fractional energy contained in 
each of these modes. For each point in Table 1, a corre- 
sponding array of energy as a function of (x, t) can be 
produced. The resolutions in x and t of each energy ar- 
ray are equal to the wavelength and period of the mode, 
respectively. Four examples are shown in Figure 4 for 
the four modes containing the largest fractional energy. 
These modes are also highlighted in boldface in Table 
1. Together these four modes comprise 64% of the total 
energy in the spectrum. Note that the left and right 
halves of the figure display properties of the westward 
and eastward propagating waves, respectively. The con- 
tours describe how energy for the resolved mode with 
wavelength and period (A X ,T) is distributed in (x,i). 
When these figures are compared to Figure lb, it can 
be seen that the wavelet analysis describes how energy 
in a given spectral mode is localized in (x,t). 


Table 1. Percent Energy in the Wavelet Spectrum as a Function of Horizontal 
Wavelength and Period 


T, min 





J 

u 

*-c 




768 

384 

192 

96 

48 

24 

12 

6 

3 

256 

0.11 

<0.01 

0.01 

0.03 

0.04 

0.06 

0.02 

<0.01 

<0.01 

128 

0.04 

<0.01 

0.02 

0.05 

0.06 

0.08 

0.03 

<0.01 

<0.01 

64 

0.04 

0.05 

0.20 

1.50 

1.73 

0.44 

0.07 

0.01 

<0.01 

32 

0.02 

0.04 

0.20 

1.46 

9.46 

6.01 

1.33 

0.08 

0.04 

16 

<0.01 

<0.01 

0,02 

0.17 

4.41 

15.83 

10.98 

0.33 

<0.01 

8 

<0.01 

<0.01 

<0.01 

0.02 

0.24 

5.52 

29.16 

5.47 

0.05 

4 

<0.01 

<0.01 

<0.01 

0.02 

0.07 

0.20 

1.25 

2.77 

0.12 
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Figure 4. Energy distributions in (x, t) from the wavelet analysis for the four modes in boldface 
in Table 1. These four modes contain the largest fractional energies in the wavelet spectrum 
(Table I). Maps like these show how spectral energy at a given (Jfe.w) is distributed in (z,t). 
Comparison to Figure lb shows how the wavelet analysis captures the spatial/temporal variability 
in the spectral modes. The contours label energy distribution as percent of maximum at 10% 
intervals . 


Spatially, Figure 4 shows that energy tends to be 
tightly constrained within narrow regions of x, The 
span of the distribution in x can also be seen to roughly 
vary with the wavelength of the mode. The relationship 
between the energy distribution in time and the wave 
period is not so obvious. To develop an objective es- 
timate of wave packet widths in space and time, the 
following method was applied to each energy distribu- 
tion in Table 1 (like those shown in Figure 4). (Note, 
however, that the energy arrays representing the longest 
resolved period and wavelength (top row and first col- 
umn of Table 1) are omitted because these arrays only 
resolve a single point in t and x } respectively, and con- 
tain very little energy anyway.) East and west halves 
of the domain are analyzed separately to separate east- 
ward and westward propagating modes. 

To estimate equivalent packet width in x , the energy 
density for each mode Ekui{x^t) is summed over time: 

^ku{ x i) — ^ ^ Ektj( x i, tj ) 

i 


where i and j are dummy indices representing the dis- 
crete grid points in the wavelet energy arrays. The wave 
packet width xwp is then estimated as the width of a 
box whose length is equal to the maximum energy in 
Ekui x i) an d whose area is the total energy in the mode, 
or 

_ St 

XWP ~ -maxlEiJxt)] 

Figure 5 illustrates this definition of wave packet width. 
The equivalent width xwp is plotted versus wavelength 
in Figure 6a. An energy weighted least squares straight 
line fit is computed and overplotted as the solid line 
in Figure 6a. The correlation coefficient is 0.91, with 
slope = 1.24 ± 0.06 and intercept = 27 ± 2 km. Fits 
to eastward and westward propagating modes were not 
significantly different, so both are included in Figure 
6. The dashed lines show the 95% confidence limits to 
the weighted linear fit. The dotted line represents the 
resolution of the wavelet analysis. These packet widths 
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Figure 5* Illustration of the definition of wave packet 
width. The vertical axis is energy per unit distance, 
and the horizontal axis is distance. Packet widths xwp 
are determined as the width of a box (dotted line) 
whose height is the maximum of the energy distribu- 
tion and whose area is the total energy in the mode 
(jb,u;). This width is approximately equal to the full 
width at half maximum of a Gaussian shaped packet 
with the same amplitude (dashed line). The cosine en- 
velope with wavelength equal to i xxwp (solid line) also 
encloses the same area and is used to reconstruct the 
wind and temperature fields at a given height as in Fig- 
ures 7 and 17. 

are fairly narrow compared to the wavelength of the 
wave they enclose, casting some doubt on the use of 
the simple oscillatory form assumed in x to compute 
derivatives in the linear equations, i.e., d/dx = ik. In 
section 4 the linear analysis will be tested against the 
full nonlinear model below 32 km where they overlap, 
and the linear treatment in fact produces quite realistic 
results. 

An exactly analogous procedure was performed to de- 
termine wave packet widths in time, twPi that are plot- 
ted in Figure 6b along with the weighted least squares 
straight line fit and 95% confidence limits. The cor- 
relation coefficient is again 0.91 with slope equal to 
1.37 zb 0.06 and intercept equal to 0.50 ± 0.02 hours. 


Modes with periods beyond 2 hours are weighted quite 
low and do not affect the fit, but there is indication of an 
upper bound to the wave packet width of ^2.5-3 hours 
for periods longer than an hour. An upper limit of 3 
hours for the wave packet width in time is then chosen 
as a modification to the straight line fit in Figure 6b. 

These results are used to define the amplitude cor- 


rection factor 

in (3): 


n* 

: — w p . jcpy p — 1.24A* + 27 km 

(5) 

N x “ 

N x Ax 

n t 

twp twp = 1.37 T + 0.50 hours 

(6) 

N t ~ 

N t At 


where A x is wavelength in kilometers, and T is wave 
period in hours. (Note also that these correction fac- 
tors cannot exceed unity, a value which would indicate 
a wave packet width equal to the size of the domain.) 
Table 2 lists values of the correction factor N x N t /n x n t 
as a function of wavelength and period. Modes con- 
taining less than .01% of the spectral energy have been 
omitted. With these correction factors and equation (3) 
the maximum vertical velocity amplitude at z =13 km 
for any single mode predicted from the power spectrum 
(Figure 2) is 1.8 ms" 1 . These amplitude corrections are 
applied for the purpose of determining breaking levels 
and wave saturation. 

3. Linear Wave Propagation With 
Saturation 

The following analysis is designed to study the ver- 
tical propagation of waves produced in the nonlinear 
squall line model to altitudes much greater than the 
current top at 32 km, into the upper stratosphere, meso- 
sphere, and lower thermosphere. The difficulties in ex- 




Figure 6. Equivalent wave packet widths in (a) x and (b) £, derived from the wavelet analysis 
plotted as a function of wavelength and period respectively. The result for each mode in Table 
1 is represented by a plus symbol. Solid lines show the energy-weighted least squares straight 
line fits to the wave packet results, which have high correlation coefficients equal to 0.91. (Note 
from Table 1 that many of the modes contain very little energy.) Dashed lines show the 95% 
confidence limits on these fits, and the dotted line in each panel indicates the resolution limit 
of the wavelet analysis. The fits define simple relationships between wave packet widths and 
frequency and wavenumber of the wave which will be used in subsequent analysis. 
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Table 2. Wave Packet Correction Factor ( N a Nt)/(n 9 nt ) as a Function 
of Horizontal Wavelength and Period 


T, min 





} 

km 




768 

384 

192 

96 

48 

24 

12 

6 

3 

256 

1.0 


1.9 

3.5 

6.0 

9.1 

12.3 



128 

1.0 


3.6 

6.6 

11.1 

16.9 

22.9 



64 

1.7 

3.3 

6.3 

11.4 

19.3 

29.4 

39.9 

48.5 


32 

2.7 

5.3 

10.0 

18.2 

30.8 

46.9 

63.6 

77.3 

86.7 

16 



14.3 

25.9 

43.8 

66.7 

90.4 

110.0 


8 




32.9 

55.5 

84.6 

114.6 

139.4 

156.2 

4 




37.9 

64.0 

97.8 

132.3 

160.9 

180.4 


tending the mesoscale simulation to study wave-mean 
flow interactions at these altitudes were described in the 
introduction. The linear propagation analysis derived 
here can be used to estimate the effects of the convec- 
tively generated waves on the mean flow for realistic 
middle atmosphere winds, and it also provides a means 
of comparing the nonlinear simulation results to recent 
observations at altitudes near the mesopause [Swenson 
and Mende , 1994; Taylor et a/., 1995]. 

Each point in the power spectrum, P(k 1 u ;), in Figure 
2 is assumed to describe the amplitude of a monochro- 
matic wave A w (k } cv) y at 13 km altitude, and centered 
at x = 420 km at the center of the storm. Note that 
the frequencies in Figure 2 are intrinsic frequencies at 
the lower boundary of the linear analysis {z — 13 km), 
where the mean wind speed is set to 16 m s -1 to match 
conditions in the nonlinear model from which Figure 2 
was derived. The power to amplitude conversion is de- 
scribed by equations (3), (4), (5) and (6). Linear grav- 
ity wave theory and the polarization relations [Gossard 
and Hooke f 1975, pp. 97-100] relate this vertical veloc- 
ity amplitude to the toted wave energy per unit density 
E [Fritts and VanZandt y 1993], 

N 2 

E=^Al, (7) 

where N is the buoyancy frequency and w is the intrin- 
sic frequency. Similarly, the square of the horizontal 
velocity amplitude A 2 can also be related to E and A ^ , 

(8) 

Each wave will follow a ray path, ( X(z )), described by 
Lighthill [1978, sec. 4.6] as 

* w =/(£),./ 2 ' « 

o 

the slope along the ray is defined by the direction of 
energy propagation relative to a stationary reference 
frame, 


/dx\ _ U(z) + dw/dk 

\ dz )r ay du) / dm 


where U(z) is the mean wind and (fc, m) is the wavenum- 
ber vector. The linear nonhydrostatic dispersion rela- 
tion can be written as 


±Nk 

(k 2 + m 2 ) 1 / 2 


= ±iVcos 9, 


( 11 ) 


where 6 is the angle at which lines of constant phase lie 
from the vertical. Equation (11) assumes that cv 2 >> 
/ 2 and m 2 >> [2H)~ 2 (/ is the Coriolis parameter; H 
is density scale height). The effects of rotation could 
ecome important when the waves are very near their 
critical levels. Marks and Eckermann [1995] point out 
the importance of the density scale height term on the 
prediction of turning point altitudes; however, for the 
spectrum considered in this work, neither of these ef- 
fects is important in estimating the mean flow forcing 
(see section 5). Differentiation of (11) with respect to 
k and m, and substitution into (10), yields 

kU 2 „ 

= =F-rrsec 2 0 esc 0 - tan 0. (12) 

ray 

In general, U, N y and 6 all vary with height. This ex- 
pression describes the slope of the ray path as a function 
of height and with a simple numerical integration using 
the trapezoidal rule will define the ray path of the wave. 
Note that the sign convention used here is iv > 0 (9 < 0) 
for eastward intrinsic phase speeds, cv < 0 (6 > 0) for 
waves with westward intrinsic phase speeds, and k pos- 
itive definite, which is different from that of Lighthill 
and accounts for the sign differences in this derivation. 
The upper sign, then, in these equations refers to east- 
ward and the lower to westward propagating waves in 
the frame of the storm. 

Two conditions on the intrinsic frequency define the 
maximum height to which waves can propagate: 

1. Where the wave phase speed equals the mean wind 
speed the intrinsic frequency cv — ► 0, 6 — > 90°, and the 
wave will be absorbed there. This level is traditionally 
referred to as the critical level. 
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2. Where the intrinsic frequency, w — ► N, and 
cos(0) — ► 1, the vertical wavenumber, m — > 0. So- 
lutions above this level have m imaginary, and so the 
wave undergoes total internal reflection. In the case 
of nonzero curvature in the background wind U(z) the 
solutions can theoretically be modified to permit prop- 
agation through this level if the condition u j = N occurs 
very near an inflection point in the mean wind profile. 
This situation is rare enough to neglect, so it will be 
treated as a simple reflection and such levels will be 
referred to as turning points. 

The ray path integration for each wave is taken only 
up to the level just below the first critical level or turn- 
ing point. At and above that height the wave energy is 
set to zero. Below the critical level or turning point, and 
in the absence of dissipation, the wave amplitude will 
change with height according to conservation of wave 
action flux along a ray [Lighthill, 1978, p. 331]: 


P Ec g* 

oj 


= const, 


where p is density, and c 3m = dw/dm is the vertical 
component of the group velocity. With (7) this can be 
rewritten in terms of the vertical velocity amplitude to 
describe the change in with height, 



( N 2 /u> 2 ) 



(13) 


The zero subscript refers to the lower boundary value. 
At some altitude this predicted amplitude will grow to 
the point where the wave would become unstable. Here 
the convective instability criterion is used to determine 
this wave breaking level, and the wave amplitude above 
that height is maintained at or below this instability 
limit. This is the traditional saturation theory, reviewed 
by Friits [1984]. At the breaking level the total Eulerian 
parcel velocity exceeds the ground relative phase speed, 
and linear theory breaks down. Dewan and Good [1986] 
derive the instability criterion for the vertical velocity, 


though this situation never occurs in the cases consid- 
ered here. 

The saturation condition (14) implies a certain shape 
of the power spectrum as a function of vertical wavenum- 
ber, m. VanZandt [1982] and many other authors since 
have described the universality of the shape of the corre- 
sponding horizontal velocity power spectrum P u {tti) in 
observations, which follows a power law approximately 
proportional to m -3 . For the vertical velocity the ana- 
log to equation (2) for the vertical wavenumber power 
spectrum and (14) imply the saturated power spectral 
density, 

P r = “ 2 

(Am) 2m 2 Am 

where Am is the spectral bandwidth. The linear disper- 
sion relation (11) simplifies in the hydrostatic approxi- 
mation (m 2 >> k 2 ) to 

0 N 2 k 2 


which is a good approximation in the short vertical 
wavelength “tail” of the spectrum where saturation oc- 
curs. Substitution then gives 


, x ^ 2jfc2 

P * (ro > = 2m*Am 


(15) 


The analogous expression for the saturated horizontal 
velocity power spectrum was derived by Dewan and 
Good [1986] and, with the definition of the vertical ve- 
locity power spectrum (analogous to (2)), can be writ- 
ten as 

w = = < I6 > 


Note, however, that unlike the Dewan and Good 
[1986] analysis, no separate assumption about the band- 
width is made here to determine the shape of the sat- 
urated spectrum as a function of vertical wavenumber, 
m. The bandwidth, Am, in this case, is dictated by the 
resolution (Aik, Auj) of the power spectrum (Figure 2) 
used to define a mode in this analysis. The hydrostatic 
dispersion relation implies 


which after substitution of (11) to eliminate the vertical 
wavenumber can be written as 


A k Aw 2 

Am = T mT *ib m 


(17) 


A? < (14) 

w — Jb 2 [(JV 2 /a/ 2 ) — 1] 

The breaking level of the wave is defined as that level 
where the amplitude equals or exceeds this saturation 
limit. Above that level, wave amplitudes will be mod- 
ified as necessary by condition (14), and (13) is then 
taken to describe the linear growth in amplitude from 
one level (subscript zero) to the next. Waves that are re- 
flected before reaching the breaking level are discarded 
and will have no net effect on the mean state. If satura- 
tion were to occur below the turning point, the wave 
would contribute to the mean state interaction only 
between the breaking level and the turning point, al- 


with A k and A lj fixed constants in this case. There- 
fore the Dewan and Good [1986] assumption, Am oc m, 
will only apply in those regions of the spectrum where 
the first term on the right-hand side of (17) dominates. 
More generally, the saturated spectra in these calcula- 
tions will follow P w (m) oc m” 5 to m -6 and P u {m) oc 
m -3 to m -4 . Numerous spectral observations of hor- 
izontal wind fluctuations show that a P u (m) oc m 3 
spectral shape is ubiquitous in the middle atmosphere 
[ VanZandt , 1982; Smith et al ., 1987]. This observed 
power law relationship and the associated amplitude 
limits have been attributed to the process of satura- 
tion. The spectra produced in this model are roughly 
consistent with observations but do not provide any new 
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insights into the origin of the observed spectral features. 
An example of the saturated P u (m) spectrum produced 
in this model will be shown in section 5. 

3.1. Energy, Momentum, and Wave Action 
Fluxes 


3.2. Wave Interactions With the Mean State 

The fluxes derived above can be used to infer the 
effects of wave dissipation on the mean state. A vertical 
gradient in the momentum flux (19) implies a drag force 
per unit mass on the mean flow U : 


Wave amplitudes derived from equations (13) and 
(14) can be used to compute the fluxes of energy and 
momentum associated with the convectively generated 
waves. These relationships have been derived by Fritts 
and VanZandt [1993] and reproduced here in slightly 
different form, in terms of the information contained in 
the power spectrum (Figure 2). By using equations (7) 
and (11), energy flux per unit density can be written as 

F e = c 9 .E = ±-(-- l) (18) 

Momentum flux per unit density is computed as the al- 
and t-averaged product of horizontal and vertical veloc- 
ity perturbations, uw, and with (8) becomes 



These expression differ from those of Fritts and Van- 
Zandt in the neglect of rotation effects (/ = 0) and in 
that no “anisotropy factor” appears here because (18) 
and (19) describe a single wave mode rather than a full 
wave spectrum. Also unlike the work of Fritts and Van- 
Zandt, no assumptions about the shape of the spectrum 
are made to integrate these expressions over all frequen- 
cies and wavenumbers. Instead the integration is per- 
formed by summing over the individual points in the 
spectrum. This method allows the spectrum from our 
peculiar storm source to be free to take any shape. The 
Fritts and VanZandt model is instead intended to cap- 
ture the more general properties of middle atmosphere 
waves. 

An expression for wave action flux in terms of the 
same variables can also be derived using equations (7) 
and (11), 



Equations (18), (19), and (20) describe fluxes associ- 
ated with a single wave mode defined by (fc,u;). Before 
computing net effects of the full spectrum of waves, the 
amplitude adjustments previously applied in (3) for the 
purpose of predicting breaking levels and saturation ef- 
fects must first be removed. This amounts to dividing 
the amplitudes by the correction factors ( N x N t /n x n t ), 
or 



The same factor must also be applied to the energy 
and wave action fluxes before any net effects can be 
computed. 


D - = ~wfi F ^ < 22 > 

if the momentum flux is averaged over a suitably large 
horizonal domain. ( p is the basic state density profile.) 
Integration over the spectrum 


Yl( D *)ku 


gives the net zonal drag on the mean flow, averaged over 
a distance of 384 km and a time interval of 4 hours, cor- 
responding to the averaging over the full (aj,t) domain 
of the original spectrum. 

Durran [1995] underscores the difficulty in using (22) 
to describe the momentum budget in numerical models 
with a limited horizontal domain. In particular, contri- 
butions to the full nonlinear momentum budget due to 
differences in the upstream and downstream values of 
(pu 2 +p) are significant. The net mean flow acceleration 
associated with breaking mountain waves in his model is 
spread over very large distances away from the region 
containing the largest wave perturbations. The same 
difficulties arise in evaluating the momentum budget in 
the nonlinear storm model. 

With this in mind, the mean flow effects derived from 
(22) in the linear analysis described here are further nor- 
malized to represent an average over the circumference 
of the earth at a latitude of 40° and over 1 day of time. 
The result given in ms" 1 day -1 represents the zonally 
averaged drag exerted on the mean flow, and averaged 
over a day’s time, due to one storm active for a period 
of 4 hours. This normalization also allows more direct 
comparison of these results to global parameterizations 
of gravity wave drag. 

Gravity wave dissipation also leads to transfer of en- 
ergy to the mean flow. The energy dissipation rate can 
be derived from conservation of wave action flux. By 
defining the wave action per unit density A = E/w and 
the associated flux of A, F \ = A c 9m , then 

|(M) = |(PF A ). (23) 


Note that both sides of (23) are zero in the absence 
of wave dissipation. Solving for the energy dissipation 
rate, e, 


€ = — 


dE 

~di 


— w d f N 


(24) 


describes the rate of energy transfer from the waves to 
the mean flow along a ray and with equation (20) can be 
computed readily from the linear analysis. Analogous 
to the net zonal drag calculation, contributions from in- 
dividual modes are summed to compute the net energy 
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dissipation rate and normalized as described above. So 
the results then given in m 2 s~ 3 represent the zonally 
averaged rate of energy dissipation to the mean state, 
averaged over a day’s time, due to waves produced by 
a single storm active over a 4-hour period. 


4. Comparison of Nonlinear and Linear 
Models 

As a test of the validity of the method and assump- 
tions within the analysis described in sections 2 and 
3, comparisons between these linear wave propagation 
calculations and the full nonlinear model between 13- 
and 32-km altitude are presented in this section. Ver- 
tical velocity uj(x), momentum flux profiles, and the 
power spectrum versus vertical wavenumber P w (m) are 
each computed from the linear model for comparison to 
the full nonlinear results. Ray paths and breaking level 
predictions from the linear model are also examined for 
qualitative agreement with the nonlinear model results. 

4.1. Reconstruction of the Vertical Velocity u;(x) 

By using only the power spectrum (Figure 2) and the 
results of the linear wave propagation analysis, realiza- 
tions of the field tu(a;) at selected heights can be recon- 
structed and compared to instantaneous distributions 
of tu(x) in the nonlinear model. These reconstructions 
represent properties averaged in time. The comparison 
provides a test of the derived wave packet widths in x , 
as well as the overall amplitudes implied by the wave 
packet width adjustments in x. 

The ray tracing analysis provides the location of each 
point in the spectrum X{z y kiyWj) via equations (9) and 
(12). Ray paths in a stratosphere with constant 16 m 
s -1 winds, like those in the nonlinear simulation, are 
straight lines and form a simple fan pattern emanat- 
ing from the central source region. The fan pattern 
appearance is similar to the pattern observable in the 
nonlinear model, but the time dependence of the wave 
forcing, especially evident in the eastward propagating 
waves (see Figure lb), is not captured in this analysis. 
The ray location at a specified altitude is taken as the 
central locus of a packet of wave energy with amplitude 
determined by wave action flux conservation/saturation 
(equations (13) and (14)). The energy in the packet is 
then distributed in x within a half-cosine envelope with 
wavelength i vxwP' This envelope was illustrated in Fig- 
ure 5 and has the same amplitude and total energy of 
the square wave packet with width x\yp. Each mode 
is then renormalized by multiplication of the tempo- 
ral width factor given by (6), returning the power back 
to time-averaged values. The bias factor used to es- 
timate the true power at the peaks in the spectrum 
( N x N t /n x n t ) is also removed before the contributions 
from each spectral mode are summed to produce pro- 
files of io(x) shown in Figure 7 for an altitude of 30 
km. Since the sum of the contributions from each mode 
may depend strongly on the phase of the wave within 
each wave packet, the phase is randomized before the 
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Figure 7. Reconstructions of the vertical velocity dis- 
tribution u>(x) at z=30 km from the spectral input and 
linear wave propagation analysis, (a)-(d) Different real- 
izations of the field created by selecting different phase 
randomizations for each wave mode contributing to the 
total. The mean flow 17= Oms -1 in these calculations 
so that results are shown in the storm-relative frame for 
comparison to the nonlinear model (Figure 8). 


sum is computed. Four different realizations of the field 
w(x) are presented in Figures 7a-7d, each computed 
with a different random number sequence to determine 
the phase randomization. Figures 8a-8d show four real- 
izations of the nonlinear model uj(x) at 30 km. Figures 
7 and 8 show similar amplitudes and similar variations 
as a function of x. The similarities suggest that the 
wave packet width in x is approximately correct and 
that the computed ray paths are capturing the trajec- 
tories of the wave energy packets fairly well (despite the 
many simplifying assumptions like the point source ap- 
proximation at 13 km and the time-averaging effect of 
using the power spectrum as input). Similar compar- 
isons at 20 km (not shown) demonstrate that the growth 
of wave amplitude with height is being accurately sim- 
ulated by the linear calculation; however, the distribu- 
tions in x appear less similar to the nonlinear model at 
this altitude much closer to the prescribed source (at 
13 km), presumably because of the inaccuracy of the 
point-source approximation. 

4.2. Momentum Flux Profiles 

Time- and as-averaged profiles of momentum flux can 
also be computed from the linear wave propagation 
analysis and compared to the nonlinear simulation. The 
nonlinear results are averaged over 4 hours and ±400 
km from the central source region, with east and west 
halves of the domain averaged separately. The results 
are shown in Figure 9a. Momentum flux in the lin- 
ear model can be calculated from equation (19), and 
contributions from each mode summed after removal of 
the wave packet and bias adjustments (21). Where ray 
paths extend beyond ±400 km from the source, those 
waves do not contribute to the average for the purpose 
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Figure 8. Vertical velocity distribution w(x) at z=30 
km from the nonlinear convection simulation, (a)-(d) 
Four times (3.03, 3.53, 6.03, and 4.87 hours) selected 
randomly for comparison to the four reconstructions 
in Figure 7. The comparison shows very similar am- 
plitudes, as well as similarities in the distribution of 
horizontal wavelength in x. This comparison provides 
a check on the computed ray paths and on the wave 
packet widths in the linear analysis. 


of comparison to the nonlinear model. Westward and 
eastward propagating modes axe averaged separately 
and overplotted in Figure 9a as heavy lines. The lin- 
ear and nonlinear results are very similar at z > 15 
km, suggesting the momentum in the vertically propa- 
gating waves is well represented by the linear spectral 
model. At lower altitudes the nonlinear model shows 
an abrupt decrease in magnitude from the tropopause 
to 15 km in the westward average. This may be an 
indication of wave dissipation in the nonlinear model 
at these altitudes, where a thin layer of shear develops 
that is not considered in the linear comparison. The 
very low frequency waves not resolved in the spectrum 
may also be contributing to the domain-averaged mo- 
mentum flux at these lower altitudes in the nonlinear 
model. Nonetheless, Figure 9a suggests the flux to the 
upper stratosphere and mesosphere is well characterized 
by the linear approximation. 

A further test of the linear analysis is shown in Fig- 
ure 9b. The thin lines again show momentum fluxes in 
a nonlinear simulation identical to that of AHD, except 
for the inclusion of strong shear equal to —2 m s _1 km -1 
in the stratospheric winds above 13 km. The strato- 
spheric shear dramatically alters the vertically propa- 
gating waves but has no discernible effect on the de- 
velopment of the storm. The heavy lines in Figure 9b 
display the momentum flux profiles predicted from the 
original wave spectrum (Figure 2) and the linear wave 
propagation analysis through a background wind with 
shear identical to the nonlinear model. The linear and 
nonlinear momentum flux profiles are very similar in 
both magnitude and form. Differences are compara- 
ble to the no shear case in Figure 9a. Predicted ray 
paths also compare quite well to the nonlinear result 


(not shown). The comparison lends further credence to 
the linear wave treatment, even in this case with fairly 
strong background shear. 

4.3. Reconstruction of the Vertical Wavenumber 
Power Spectrum 

As a further test of the linear model assumptions the 
spectrum P w (m) can be computed from the spectrum 
in Figure 2, P w (k,(j), with the help of the linear grav- 
ity wave dispersion relation (11). The result can be 
compared to the spectrum computed directly from the 
vertical velocity field in the nonlinear simulation with- 
out stratospheric shear. That spectrum, described by 
AHD is reproduced in Figure 10 and compared to the 
linear model reconstruction. 

For the reconstruction the value of m at each point 
in the P(k>w) spectrum is computed from (11). The 
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Figure 9. Profiles of momentum flux pquw averaged in 
x and t. Thin lines show profiles averaged west (solid) 
and east (dashed) of storm center from the nonlinear 
convection simulation. Heavy lines show the same quan- 
tities reconstructed from the spectral input with the re- 
sults of the linear ray tracing model, (a) Comparison of 
the AHD model and the linear analysis, (b) Compari- 
son to a nonlinear simulation identical to AHD except 
with —2 ms -1 km -1 shear in the stratosphere above 13 
km. The thin lines are derived from the same spectral 
input as in Figure 9a but are now affected by the strong 
background shear. 
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Figure 10. Power spectrum of the vertical velocity 
as a function of vertical wavenumber P w (m). Thin 
lines show spectra computed directly from the nonlinear 
convection simulation (AHD), while heavy lines show 
Pw{™) reconstructed from the spectral input (Figure 2) 
by using linear gravity wave theory. This is the no-shear 
case. The solid and dashed lines represent the spectrum 
of waves west and east of storm center, respectively. 


power at each point P(k ) a>) is then binned in a verti- 
cal wavenumber array, identical to that in the original 
spectrum: 


m = 


n 

N z Az 


n- 0,1,..., 


N z 


- 1 


where N z — 128, and A z — 148 m (AHD). Power falling 
into each bin is integrated over k and lj: 


P(m)Am = P(k(m) f iv(m))AkAcv 

k t cj 

P(m) is the resulting spectrum overplotted in Figure 
10. The similarity of the result to the nonlinear model 
spectrum reflects the accuracy of applying linear gravity 
wave theory to these motions. It also suggests that the 
errors associated with truncating the power spectrum at 
Pmin are small but evident as an absence of power in the 
low-energy tail of the reconstructed P w (m) spectrum. 


5. Wave Interactions With Realistic 
Background Wind Profiles 

With wave properties obtained from the nonlinear 
convection simulation, the linear ray tracing and wave 
action conservation model can predict the behavior of 
these vertically propagating waves through climatolog- 
ical background wind and buoyancy frequency profiles. 
The CIRA [Fleming et a/., 1990] model at 40°N lati- 
tude provides background state profiles for this analy- 
sis. Profiles for April and June monthly averages are 
chosen to illustrate the nature of the wave-mean flow 
interactions. 


5.1. Breaking Levels and Critical Level Filtering 
Associated With Variable Background Winds 

The intrinsic frequency in the linear dispersion rela- 
tion (11) will vary with height when there is shear in the 
background wind as u) = too — kU(z) } so in the presence 
of shear, vertically propagating waves may find critical 
levels as described in section 3. This Doppler shifting 
also affects wave amplitude growth, breaking levels, and 
saturation of the waves via the frequency dependences 
in equations (13) and (14). 

Figures 11a and lib illustrate some of these effects 
in the linear model. The solid lines represent climato- 
logical wind profiles for April and June. Symbols mark 
the ground-relative phase speeds of the waves in the 
spectrum along the abcissa. Dot symbols mark turning 
point altitudes where waves are reflected as the magni- 
tude of their intrinsic frequency reaches the local buoy- 
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Figure 11. Zonal winds, breaking levels, and turn- 
ing points for (a) April and (b) June, with background 
states defined by the CIRA model at 40° N latitude. 
Solid lines are zonal wind profiles ( U ). Plus symbols 
mark breaking levels as a function of ground-relative 
phase speed (cc) for each wave mode. Dots denote turn- 
ing points. Critical levels are not explicitly shown but 
occur at altitudes where c* — U . Wave modes con- 
tribute to the mean flow forcing between their breaking 
and critical levels. 
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ancy frequency. Wave reflection tends to occur in re- 
gions where the background shear causes growth in the 
intrinsic wave frequency. For the climatological winds 
considered here, none of these modes saturate below 
their reflection levels, so they are discarded from fur- 
ther analysis, having no net effect on the background 
state. Plus symbols mark breaking level altitudes as 
a function of phase speed. Wave breaking tends to oc- 
cur where the background shear causes the intrinsic fre- 
quency to shrink in magnitude. If the shear is strong, 
as in the June case in the stratosphere, the breaking 
level tends to appear just below the critical level where 
w — > 0. This also tends to be true in the stratosphere 
where wave amplitudes predicted from wave action con- 
servation (13) have not grown large enough to saturate 
(14) until their intrinsic frequencies become very small. 
In the much rarer atmosphere in the mesosphere, waves 
more often break well below any critical level and so 
contribute far more energy and momentum to the mean 
state forcing via the effects of saturation. Wave criti- 
cal levels are not explicity shown but occur where the 
ground-relative phase speed matches the speed of the 
background wind. 

5.2. Saturation Effects on the Power Spectrum 
Versus Vertical Wavenumber 

The effects of saturation on the shape of the horizon- 
tal velocity perturbation power spectrum versus vertical 
wavenumber have been described by numerous authors 
[ Dewan and Good , 1986; Smith et al. } 1987; Friits and 
VanZandt , 1993] and have been offered as an explana- 
tion for the ubiquitous observations of a power curve 
with a shape proportional to m~ 3 as described in sec- 
tion 3. Fritts and VanZandt [1993] describe this satu- 
rated spectrum as 


P u (m) 


(p+ l) iV 2 

20 to 3 


(25) 


with p = |. This theory also explains the nearly con- 
stant amplitude of the spectrum with height despite the 
large changes in density from the troposphere to the 
thermosphere. In regions where waves are saturating, 
the linear wave propagation model developed here dis- 
plays characteristics similar to the theory and to spec- 
tral observations. 

The region between ~20 and 30 km in the June case 
(Figure 11b) provides an example. At these altitudes, 

waves with phase speeds between ' 20 and 0 m s” 1 

break, and their intrinsic frequencies shrink rapidly to- 
ward zero as the waves approach their critical levels. 
Alternatively, waves with phase speeds larger than 16 
m s~* are Doppler shifted to higher frequency. The 
effect of the shear on the vertical wavenumber spec- 
trum is shown in Figure 12. The thin line shows the 
original spectrum, P u (m), at 13- km altitude. Waves 
with phase speeds of <16 m s _1 are Doppler shifted 
to higher vertical wavenumbers and saturate. The east- 
ward propagating waves with phase speeds of > 16 m s -1 
are Doppler shifted to lower m and grow in amplitude 
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Figure 12. Saturation effects on the horizontal veloc- 
ity power spectrum P u (^)» The thin line shows the 
input wave spectrum at 13 km including both east- 
ward and westward propagating waves. The heavy solid 
line shows the “saturated spectrum” between 20 and 
30 km in the June case. Eastward propagating waves 
are Doppler shifted to long and grown in amplitude, 
while westward waves get shifted to short X z and satu- 
rate. The dashed line is the theoretical curve, equation 
(25). 

as they propagate vertically. The saturated spectrum 
averaged between 20 and 30 km is shown as the thick 
solid line in Figure 12, and the theoretical saturation 
limit (25) is plotted as a dashed line. The spectrum is 
noisy, since it is computed as a sum of the discrete set 
of modes, but approximately follows an m -3 to ra -4 
power law, as expected from equations (16) and (17), 
and the amplitude is similar to that given by the the- 
ory (25). (Note that variations in wave packet travel 
times associated with variations in group velocity are 
neglected in these calculations.) It should be empha- 
sized that the fit of the P u (m) spectrum in this model to 
the saturation theory (25) is a direct result of the spec- 
tral treatment (17) coupled to the imposed saturation 
condition (16). This result provides no new insight into 
the physical mechanisms responsible for the observed 
m~ 3 power law. It has been suggested that wave- wave 
interactions in fact produce the observed spectral shape 
[ Weinstock , 1990; Hines , 1991], and these have been ne- 
glected in the present model. Figure 12 shows, however, 
that the linear spectral model with saturation described 
here is at least not inconsistent with the observations. 

5.3. Zonal Drag and Dissipation Rates 

Wave drag and dissipation rates for the April and 
June cases are shown in Figures 13 and 14. The 
rates plotted represent the contribution to the mean 
flow acceleration or wave energy dissipation rate result- 
ing from the single storm of 4-hour duration averaged 
over the latitude circle at 40° N and over a day’s time. 
Also shown for reference are the parameterization re- 
sults from Fritts and Lu [1993] plotted at one-tenth 
magnitude as dashed lines. In both figures, lower panels' 
display stratospheric altitudes and upper panels show 
mesospheric altitudes. 
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Figure 13. Zonal mean drag force per unit mass on the 
zonal wind from the convectively generated wave spec- 
trum (solid lines). The dashed lines show the Fritts and 
Lu [1993] parameterization at one tenth the magnitude 
for comparison. Upper panels show mesospheric alti- 
tudes, while lower panels show the stratospheric results 
at finer scale. 


The storm forcing is not very similar in form to the 
parameterization. These high-frequency waves, forced 
in the upper troposphere, have only weak effects in the 
lower stratosphere. In the more rarified upper strato- 
sphere and mesosphere their effects can be substantial. 
The storm-generated wave effects dwindle in the lower 
thermosphere: Most of the original spectrum has been 
reflected or absorbed at critical levels below these alti- 
tudes. 

In the June case, stratospheric shear in the mean 
winds eliminates much of the spectrum carrying nega- 
tive momentum flux (waves with c x < 16 m s _1 ). These 
waves are absorbed at critical levels in the stratosphere 
(Figure lib). The critical level absorption, however, 
contributes little to the mean flow forcing (Figures 13 
and 14), because as the intrinsic frequency approaches 
zero at the critical level, so the wave energy also goes to 
zero, conserving wave action flux. Wave saturation due 
to the density effect has a far greater effect on the mean 
state. In the mesosphere, in this case, waves carry- 


ing eastward momentum flux (c x > 16 m s _1 ) saturate 
well below their critical levels and deposit substantial 
momentum and energy at these levels. The Fritts and 
Lu parameterization predicts a peak drag of ~100 m 
s -1 day -1 at 85 km in general agreement with Holton's 
[1982] original summer soltice estimate. The storm drag 
peaks at nearly 40 ms -1 day -1 and so is a substantial 
fraction of the zonal mean in the upper mesosphere. 
In April the shear is weaker, and the storm-generated 
waves create drag of opposite sign as that of Fritts and 
Lu in the mesosphere. Referring to Figure 11a, the 
waves that are breaking above 60 km have eastward 
phase speeds, so they produce an eastward drag force. 
The parameterized drag, on the other hand, has some 
tendency to oppose the mean wind regardless of such 
phase speed preferences. 

The largest uncertainty in these results lies in the de- 
termination of wave packet widths that influence where 
the waves break. To evaluate this uncertainty, wave 
packet widths defined by the 95% confidence limits on 
the straight line fits in Figure 6 are alternately sub- 
stituted for equations (5) and (6), and the wave/mean 


April 

1001 


90 



60 


50 f ... i ... i 

0.00 0.02 0.04 0.06 0.08 
Dissipation Rate (rr^s -3 ) 


April 



0.00 0.05 0.10 0.15 0.20 
Dissipation Rate (10“ 3 m 2 s~ 3 ) 


June 



0.00 0.02 0.04 0.06 0.08 
Dissipation Rate (m 2 s -3 ) 


June 

60 f — 


50 


20 > : 
io L i 

0.00 0.05 0.10 0.15 0.20 
Dissipation Rote (10~ 3 m 2 s -3 


Figure 14. Wave energy dissipation rates from the 
convectively generated wave spectrum (solid lines) are 
compared to one tenth the magnitude of the Fritts and 
Lu [1993] parameterization (dashed lines). 
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Figure 15. The effect of uncertainty in the wave energy 
packet widths on the mean flow acceleration. Solid lines 
are the same as those in Figure 14. Dashed lines show 
the effect of packet widths defined by the upper limit 
fits shown in Figure 6. Dot-dashed lines Bhow effects of 
the lower limits. 


Figure 16. The effect of uncertainty in the wave en- 
ergy packet widths on the energy dissipation rate pro- 
file. Solid lines are the same as those in Figure 15. 
Dashed lines show the effect of packet widths defined 
by the upper limit fits shown in Figure 6. Dot-dashed 
lines show effects of the lower limits. 


interaction is recomputed. The results are shown in 
Figures 15 and 16. Note that these changes do not 
result in a simple scaling of the profiles in Figures 13 
and 14. Instead the forcing is redistributed in z. The 
changes affect the breaking levels of the waves but not 
the net energy and momentum they carry at the lower 
boundary, which are fixed by the integrated power in 
the input spectrum. The changes can be significant at 
some altitudes but do not affect the overall conclusions. 
Still this uncertainty underscores the need for caution 
in the quantiative use of the forcing profiles derived in 
this analysis. 

5.4. Mesopause Perturbations 

Figure 17 shows estimates of temperature perturbs 
tions at the mesopause constructed in a manner anal- 
ogous to the vertical velocities in Figure 7. Tempera- 
ture perturbations are derived from the wind perturba- 
tions via the polarization relations [Gossard and Hooke , 
1975] and the ideal gas equation of state. Here, x = 0 
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Figure 17. Mesopause temperature perturbation 
reconstructions, £T(sc), shown as percent of the 
mesopause temperature for the (a) April and (b) June 
cases. 
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represents the location of the wave source specified at 
13 km. Filtering of the westward propagating waves 
by wind shear is particularly evident in the June case, 
where there is an absence of wave perturbations west of 
the storm source. Amplitudes to the east of the source 
are limited by saturation in this region, where positive 
shear in the zonal wind is present for both the April and 
June cases. Peak perturbations are approximately 10- 
15% in the April case, and 20-30% for June. Vertical 
wavelengths at the mesopause for these waves cluster in 
the 5- to 20-km range for April, 10- to 30-km for June. 
These are large perturbations, and they would be read- 
ily visible in observations of the OH airglow like those 
reported by Swenson and Mende [1994] and Taylor et 
al. [1995]. A pattern of increasing horizontal wave- 
length with increasing distance from the source region 
is a common feature of these calculations and is con- 
sistent with the aircraft-based observations reported by 
Swenson and Mende [1994] in the vicinity of a Pacific 
storm system. Radiative damping, not included in this 
analysis, might reduce the amplitude of the perturba- 
tions derived here. Rough calculations using the scale 
dependent damping rates of Fels [1984] suggest that am- 
plitudes in the layer between 55 and 85 km would be 
reduced by only 1% or less over most of the spectrum, 
but a few of the modes may be significantly damped 
10 %). 

6. Summary and Conclusions 

This paper describes a technique for evaluating the ef- 
fects of simulated convectively generated gravity waves 
on the middle atmosphere. The method employs linear 
gravity wave theory in a ray tracing calculation from 
the tropopause to the lower thermosphere. Conserva- 
tion of wave action flux predicts wave amplitude growth 
with height along the rays, and a saturation condition 
derived from the convective instability criterion is em- 
ployed to predict breaking levels and limits to the am- 
plitude growth. 

In the work of AHD an analysis of the gravity wave 
spectrum observed above the deep convection in their 
squall line simulation established that some of the defin- 
ing characteristics of the spectrum result directly from 
the wave-forcing mechanisms active in the model and 
that some of these characteristics have also been seen in 
ground-based observations of motions above deep con- 
vective storms. The spectrum of simulated waves above 
the storm AHD describe provides the lower boundary 
input to the ray tracing model. Input wave amplitudes 
as a function of frequency and horizontal wavenumber 
are determined by combining the results of Fourier and 
wavelet spectral analysis techniques. The Fourier anal- 
ysis provides a high-resolution two-dimensional spec- 
trum of wave amplitudes as a function of frequency 
and wavenumber, but these amplitudes are averaged 
over the entire spatial/temporal interval of the analy- 
sis. Since the wave modes observed in the convection 
simulation are not uniformly present in space and time, 
an estimate of the distribution of wave energy is pro- 


vided by the spatial and temporal information retained 
in an orthogonal wavelet analysis, although at greatly 
reduced spectral resolution. This information on the 
distribution in (x,t) of the spectral energy E(k ) w) is 
cast into an estimate of the wave energy packet widths 
in space and time and then used to adjust the am- 
plitudes in the power spectrum derived from Fourier 
analysis to represent local, rather than averaged, val- 
ues. The adjusted amplitudes allow realistic estimates 
of wave breaking levels in the ray tracing calculation; 
however, the wave packet adjustments are removed be- 
fore evaluation of any net wave/mean state interactions 
to properly conserve energy. The model is tested for 
consistency with the full nonlinear simulation below 32 
km, where they overlap. 

This ray tracing model is subsequently used to ex- 
amine the effects of realistic mean wind profiles on the 
convectively generated gravity wave spectrum. The ef- 
fects of variable buoyancy frequency are also included. 
Critical level absorption, wave reflection, and the ef- 
fects of Doppler shifting and saturation on the power 
spectrum versus vertical wavenumber are all treated in 
fairly realistic manner and produce results consistent 
with observations and theory. 

The results are cast in terms of the drag force per unit 
mass on the mean flow and wave energy dissipation rate 
and compared to the Fvitis and Lu [1993] parameteri- 
zation of gravity wave effects. The comparison suggests 
that single large storm systems could contribute signif- 
icant fractions of estimated zonal mean wave forcing at 
some altitudes, particularly in the mesosphere. Further 
theoretical work on gravity waves generated by deep 
convection will establish the generality of these conclu- 
sions. 

This model may also be useful in assessing relation- 
ships between storm generated waves and observations 
of gravity wave induced brightness variations in the high 
altitude airglow such as those obtained during the re- 
cent ALOHA 93 campaign and reported by Swenson 
and Mende [1994]. Predicted temperature fluctuations 
at mesopause altitudes derived from this theoretical 
analysis would be easily detectable in such OH airglow 
observations, and some of the characteristics of waves 
associated with the Pacific storm in the observations 
are qualitatively matched in these theoretical results. 

The model developed here also provides a means of 
comparing other convectively generated gravity wave 
simulations and may allow the assessment of the rela- 
tionship between general characteristics of convectively 
generated waves and their global effects on the middle 
atmosphere. It can also be used to test simpler parame- 
ter^ ations of convectively generated waves against cal- 
culations with more detailed spectral information such 
as the example in this work. 
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